Quantum quench dynamics and population inversion in bilayer graphene 
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The gap in bilayer graphene (BLG) can directly be controlled by a perpendicular electric field. 
By tuning the field through zero at a finite rate in neutral BLG, excited states are produced. Due to 
screening, the resulting dynamics is determined by coupled non-linear Landau-Zener models. The 
generated defect density agrees with Kibble-Zurek theory in the presence of subleading logarithmic 
corrections. After the quench, population inversion occurs for wavevectors close to the Dirac point. 
This could, at least in principle, provide a coherent source of infra-red radiation with tunable spectral 
properties (frequency and broadening). Cold atoms with quadratic band crossing exhibit the same 
dynamics. 
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I. INTRODUCTION 

Charge carriers in bilayer graphene (BLG), which con- 
sists of two atomic layers of crystalline carbon, com- 
bine non-relativistic " Schrodinger" (quadratic disper- 
sion) and relativistic "Dirac" (chiral symmetry, unusual 
Berry phase) features. Due to their peculiar nature, BLG 
holds the promise of revolutionizing electronics, since 
its band gap is directly controllable by a perpendicu- 
lar electric field over a wide range of parameters 1 - - — (up 
to 250 meV— ), unlike existing semiconductor technology. 
Moreover, unlike monolayer graphene (MLG), whose ef- 
fective model (the Dirac equation) was thoroughly stud- 
ied in QED and relativistic quantum mechanics, under- 
standing the low energy properties of BLG is a new chal- 
lenge. 

Tuning the gap through zero in BLG in a time de- 
pendent perpendicular electric field parallels closely to 
a finite rate passage through a quantum critical point 
(QCP): as the gap closes, activated behaviour and a fi- 
nite correlation length give way to metallic response and 
power-law correlations, as in a sweep through a QCP. 
During the latter, defects (excited states, vortices) are 
produced according to Kibblc-Zurck theory^. When 
the relaxation time of the system, which encodes how 
much time it needs to adjust to new thermodynamic con- 
ditions, becomes comparable to the remaining ramping 
time to the critical point, the system crosses over from 
the adiabatic to the diabatic (impulse) regime. In the 
latter regime, its state is effectively frozen, so that it 
cannot follow the time-dependence of the instantaneous 
ground states - as a result, excitations are produced^. 
Evolution restarts only after leaving the diabatic regime, 
with an initial state mimicking the frozen one. The the- 
ory, general as it is, finds application in very different 
contexts in physics, ranging from the early universe cos- 
mological evolution^ through liquid 3,4 He^ i 10 ' 11 and liq- 
uid crystal s 12 ' 13 to ultracold gases 1 ^, verified for both 
thermodynamic and quantum phase transitions 1 ^. The 



relative case of manipulating the gap - in particular in 
real time - via a spatially uniform external electric field, 
which can therefore play the role of a (time dependent) 
control parameter, establishes BLG as an ideal setting for 
the study of quantum quenches with sudden, continuous 
or any other sweep protocols 1 ^—. This in turn leads to 
the question: what might such states be useful for? 

This complex of questions is addressed here. In par- 
ticular, we compute the defect (excited state) density af- 
ter a slow, non-adiabatic gap-closing passage in BLG via 
Kibble- Zurek^ theory, taking screening between the lay- 
ers into account. The presence of excited states after such 
a quench leads to population inversion for wavevectors 
near the Dirac point in BLG (see Fig.[T]), evidenced by the 
dynamic conductivity. This could in principle provide a 
coherent source of infra-red radiation with tunable spec- 
tral properties (frequency and broadening), determined 
below in an idealised model. This is promising as there 
are only few materials that generate light in the infrared 
with tunable frequency; BLG with its unique properties 
might represent the first step towards new lasers for this 
regime. 

II. HAMILTONIAN, TOPOLOGICAL 
PROPERTIES 

We study the problem in a more general setting of 
a general class of low energy Hamiltonians, comprising 
mono- and bilayer graphene, which exhibit quantum crit- 
ical behaviour, as 

n _ ( A cj{ Px -ipy) J \ m 

\cj(p x + ip v ) J -A )' (i) 

where J is a positive integer. The energy spectrum is 
given by E±(p) = ±a/A 2 + s 2 (p) with e(p) = c,j\p\ J the 

gapless spectrum, \p\ = \Jp^ + p 2 y with spatial dimension 
d=2. 
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FIG. 1: Reversing the applied perpendicular electric field +E 
in half-filled BLG (left) at a finite rate f/r leads to excited 
states in the upper branch in accordance with the Kibble- 
Zurek theory of non-equilibrium phase transitions (right). 
The momentum distribution increases from red (0) to blue 
(1) in the spectra. Realistic quenching times provide an effec- 
tive population inversion with little effect on the layer charge 
asymmetry. 

The critical exponents can straightforwardly be read 
off. The correlation length follows from dimensional anal- 
ysis: £ - ft(c J /|A|) 1 / J , defining v = 1/J. The Hamil- 
tonian contains the Jth spatial derivative ( Jth power of 
p), which leads to z = J. The resulting scaling relation 
zv = 1 is in agreement with a linearly vanishing gap A. 
To understand the nature of this criticality, let us take 
a closer topological look at Eq. ([I) by evaluating the 
Berry curvature (Q p ) for a given J, which is related to 
the phase picked up during an adiabatic excursion in the 
Brillouin zone as^ 

n p = Vp x A(p) (2) 

with A(p) = — i(np|V p |np), |np) is the eigenfunction in 
the nth band. For Eq. ((T|) , we obtain for the z component 
of the Berry curvature per valley and spin 

A fds(p)\ 2 _ AJV>| 2 (^> 
p 2E+(p) \ d\p\ J 2(A 2 +e 2 (p)) 3 / 2 ' 1 ' 

and its integral defines a topological invariant^ as 

Cj = hj d2p ^ = i si s n ( A )- ( 4 ) 

Therefore, the sign change of A corresponds to a change 
in the topological properties of Eq. (TTJ) . In addition, the 
Hall conductivity also exhibits a step as A passes through 
zero, and depends on the very same topological invariant 
per spin and valley as 

e 2 e 2 J 
a xy = —Cj = —-sign(A). (5) 

States with different values of Cj can be regarded as be- 
longing to distinct phases, similarly to the a xy plateau 



phases of the integer quantum Hall effect^. Note, that 
low energy Hamiltonians like Eq. ((T|) usually occur pair- 
wise (i.e. at the K and K' points in the Brillouin zone 
for graphene). Therefore the topological invariants, Cj 
from different valleys, add up to integer (not necessar- 
ily zero) Chern numbers. A A from spin-orbit coupling 
can trigger a non-zero Chern number, while the contribu- 
tion from different valleys due to a staggered sublatticc 
potential or bias voltage lead to a zero Chern number, 
although each valleys can have non-trivial topology with 
finite Cj. Spin degeneracy also leads to an additional 
factor of 2. 



III. QUENCHING THE GAP 

We are interested in the quantum quench dynamics 
when the gap varies as A(t) = Arji/r (up to logarithmic 
corrections, as analyzed below) and t £ [— r, r]. Accord- 
ing to Kibblc-Zurck scaling^, the resulting defect (extra 
electron/hole on the hole/electron rich layer, respectively, 
equivalent to excited states in the upper branch in this 
case22) density is p ~ T -dv/{xv+i) ^ wn j cn leads to 

P ~(A /t) 1/j . (6) 

The matrix structure of Eq. (TTJ) allows us to connect 
our problem to the Landau-Zencr (LZ) dynamics 2 ^ by 
analysing the solution of 

ihd t ^{t) = H^(t), #(-t) = (7) 

where H*f>± = E±fy±, and the quantity of interest 
is ^{t). Considering finite temperatures amounts to 
change the initial condition as a combination of posi- 
tive and negative energy states. However, as long as 
fcgT Ao, our results hold. The exact solution for 
the diabatic transition probability between final ground 
and excited states at momentum p for e{p) <C Ao gives 
for the momentum distribution of excited states in the 
upper branch (Fig. [T]) and and the resulting total defect 
density 

P p =cxp(-7re 2 (p)r/^A ) , (8) 
Ac f ,2 p A C T(1/J) (HA \ 1/J 

per valley, spin and unit cell, with A c the unit cell area. 
This agrees with Kibble-Zurek scaling in Eq. (J5|). How- 
ever, the present approach also provides the explicit nu- 
merical prefactor for arbitrary J, similarly to the quan- 
tum Ising model2i Note that the bigger J, the larger 
(and the more insensitive to r) the resulting defect den- 
sity, on account of the larger the number of low energy 
states (uj 2 / j ) within an energy window u> around the 
Dirac point. 

Since the number of defects from Eq. ^ progressively 
increases with decreasing r, it is important to address its 
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validity. From Ref. the borderline between a sudden 
and slow quench is determined from hdA/dt ~ A 2 , which 
yields tAq ~ h. Thus, our results apply in the slow 
quench regime when tAq > h, while the sudden quench 
region sets in for rA < h. 



IV. PHYSICAL REALIZATION 

A. Monolayer graphene 

The J = 1 case with c\ = vp « 10 6 m/s is realized in 
MLG — , where the spinor structure encodes the two sub- 
lattices of the honeycomb lattice. The control or even 
the very existence of a gap there remains an open is- 
sue. Dirac fermions with linear band-crossing can alter- 
natively be realized in optical lattices^ 6 -, where the on-site 
energies of different sublattices are under control, allow- 
ing for the introduction of a time dependent mass gap. 
The quantum-Hall step of Eq. ((SJ) represents the hall- 
mark of a single Dirac cone. 



B. Bilayer graphene with screening 

The J = 2 case with C2 = l/2m (m ps 0.03m e ) co- 
incides with the low energy Hamiltonian of BLG— for 
energies below tj_/4, with t± ~ 0.3 — 0.4 eV the inter- 
layer hopping, and the spinor springs from the two lay- 
ers. Keeping BLG at charge neutrality by cither isolating 
it from the rest of the world in a perpendicular electric 
field, or by using a dual-gate structure^— >2£, a continuous 
change of the gate voltage results in closing and reopen- 
ing the gap, as the density imbalance between the layers 
is inverted. Additionally, it also changes the topological 
properties of the model as in Eq. ([3]). However, screening 
due to electron interactions becomes relevant in this case, 
and the induced gap is related to the external potential, 
U e xt ao 2 ' 29 



2A = U e 



e 2 dSn 
2A c e r so ' 



(10) 



where Sn = J2 P ( n ip ~ n ^p) IS the dimensionless density 
imbalance between the two layers with n,i P the particle 
density of state p on the ith layer. In equilibrium, to a 
good approximation, the induced gap is given byii 2 - 



A 



1 + A In 



4tx 

\Uext 



(11) 



and the density imbalance reads 

5n = 4/> Alii(|A|/2tj.) J (12) 

with A = e 2 dpo/A c e r eo ~ 0.1 — 0.5 the dimensionless 
screening strength, d w 3.3 A the interlayer distance, £o 
the permittivity of free space and po = A c m/2Trh 2 the 
density of states per valley and spin in the limit A — > 0. 
For Si02/air interface, e r « 2.5 (e r = 25 for NH3, e r = 
80 for H2O), which reduces the effects of screening. 



V. NUMERICS 

In a quench of a time dependent external potential 
in BLG, the induced gap couples the two- level systems 
(stemming from the 2x2 structure of Eq. (fTJ, labeled 
by p) via the Sn term in Eq. (|10l) . The problem would 
require the solution of a continuum of coupled differential 
equations, which is not easy, even approximately. We 
mention that the case of a single level (only one p mode), 
in which case Sn = n\ v — n% v in Eq. (|10p . is known as 
the non-linear LZ model^ ., and the resulting dynamics 
differ qualitatively from the conventional one, possessing 
nonzero transition probability even in the adiabatic limit 
for strong non-linear coupling. 

The analysis is simplified considerably by the observa- 
tion that a single level cannot have a strong impact on the 
dynamics of the others due to the large number of terms 
in the sum for Sn. Thus, it looks natural to replace the 
non-linear term by an average density imbalance, inde- 
pendent of the explicit time dependence of n\ p (t) — n2 P (t) 
for a given p, hence decoupling the LZ Hamiltonians for 
distinct p's. 

When U ext changes fully adiabatically, the result- 
ing gap and density imbalance are given by Eqs. (|11[) 
and (|T2"j) , respectively. For slow, nearly adiabatic tempo- 
ral changes of the potential, only a small fraction of terms 
in the Sn sum is expected to behave truly diabatically 
(contribution from states nearest to the gap edges). Thus 
we assume that the gap is still given by Eq. (fTTj). and 
establish self-consistency by verifying that the resulting 
density imbalance satisfies Eq. (TT2"|) . Although the usage 
of Eq. (TTT)) simplifies the picture, it still differs from the 
conventional LZ form, i.e. subleading logarithmic terms 
are inevitably present albeit with a reasonably small pref- 
actor A. Fortunately, one can invoke the extension of the 
Kibble-Zurek mechanism for non-linear quenches to esti- 
mate the resulting defect densit y 16 ' 17 (note the difference 
between a non-linear quench on the LZ proble m 16 ' 17 and 
the non- linear LZ problem^). The logarithmic terms in 
Eq. pT|) can be considered as "zeroth" powers, therefore 
the resulting quench is still " linear" , with subleading log- 
arithmic corrections. 

The inset of Fig. [5] shows the density imbalance, ob- 
tained from solving numerically the LZ problem (Eq. ([7])) 
with the adiabatic screening potential (Eq. (|11[0 for BLG 
with a linearly varying external potential, 



Uext(t) = U t/T, 



t G -t, r 



(13) 



The numerical results arc compared to those from 
Eqs. (fTTj) and (TT2")k the imbalance is rather well described 
by the equilibrium, fully adiabatic (r — > 00) expression 
(dashed- green line), therefore our decoupling of the cou- 
pled non-linear LZ problem by the adiabatic potential 
for slow enough quenches with Eq. (fTTj) works satisfac- 
torily. This validates our average field decoupling proce- 
dure. Note that to the density imbalance in Eq. (|12l) all 
states up to the cutoff, t±, are contributing. On the other 
hand, defect production occurs at very low energies, close 
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to the touching point of the gapless branches, whose con- 
tribution to the imbalance is negligible in the limit of the 
size of the initial gap, Eq. (TIT]), A A = A|[/ Mt =[/ < tj_. 

The number of defects (excited states in the upper 
branch) created in an external potential, U ex t(t) = 
Uot/r, t e [— t, r], follows Eq. © even in the presence 
of screening as 



(14) 




tA ' 



where A = |[/ /2|, A A = A\ Ueait=Uo . Eq. dHJ) to- 
gether with Eq. (|9|) are the central results of our Kibblc- 
Zurek analysis. The numerical data fitted with p/ pqAq = 
C(^-) Q /2, and both the prefactor C and the exponent a 



are compared to the expected values, namely y A A /Ao 
for the coefficient and 1/2 for the t exponent for vari- 
ous values of A, summarized in Table U and shown in 
Fig. [2J The agreement is indeed remarkable, the slight 
mismatch in the exponent 1/2 being due to the sublead- 
ing logarithmic terms in Eq. (| 1 1|) for stronger screening. 
Since A po ~ 1CT 3 for A ~ ij_/10, the resulting density 
of defects per unit area (including spin and valley) falls 
into the order of y/h/rA n x 10 12 cm -2 , and can take the 
value 3 x 10 9 cm -2 for quenching time r ~ 1 ns, cor- 
responding to a ramping rate Ao/r ~ 10 7 cV/s. Note 
that this density corresponds to the electrons/holes in 
the otherwise empty/occupied upper/lower branch, and 
does not by itself imply any particular real space den- 
sity modulation, since these states contribute negligibly 
to the layer charge imbalance. A moderately slow quench 
implies tAq/K ~ 10 — 100 with Ao ~ t±/10, translating 
to t ~ 0.1 — 1 ps. Different non-linear sweep proto- 
cols^— lead to similar conclusion: the steeper (more 
non-adiabatic) the quench, the bigger the defect density 
produced. 

Our results are robust with respect to variations in 
the band structure, e.g. extra hopping terms or large 
asymmetry gap. The quadratic spectrum of BLG with 
J = 2 changes to linear one (J = 1) at the vicinity of 
the Dirac point (~ 10 K range) due to trigonal warp- 
ing, which could affect the scaling of the defect density 
(1/s/t — >• 1/t) for slow quenches. Excitonic effects can 
either renormalize the gap in biased BLG or open small 
gaps in unbiased BLG^I, which can be overcome by the 
electric field without affecting our findings. 



A 





0.1 


0.2 


0.5 


v/Ax/Ao from Eq. 


1.00 


0.88 


0.80 


0.64 


a/Aa/Ao from the fit 


1.00 


0.87 


0.78 


0.64 


exponent (a) 


0.50 


0.52 


0.53 


0.55 



TABLE I: The numerically obtained values of the coefficient, 
^/Aa/Ao and the exponent 1/2 of the defect density from 
Fig. [2]for t±_ — 5Uo, compared to the values based on Kibble- 
Zurek scaling and Eq. (|lip . 




FIG. 2: (Color online) The density of defects created dur- 
ing the quench per spin, valley and unit cell in BLG with 
screening is shown for U ex t = Uot/r, t± — 5Uo, A = (blue, 
circle), 0.1 (red, square), 0.2 (black, triangle) and 0.5 (green, 
star) from top to bottom. The symbols denote the numerical 
data, the solid lines are fits using p/poAo = v(-x~) Q '• The 
inset shows the time dependent density imbalance of BLG 
per spin and valley in a linear external potential with strong 
screening (A = 0.5) with tAq/H = 1 (blue), 10 (red) and 100 
(black) from top to bottom. The green dashed line shows the 
fully adiabatic (equilibrium) result with r — >■ oo, Eqs. (|ll|l - 
(|12p . which is approached fast with increasing r. Given the 
simplicity of our self-consistent average field procedure, the 
agreement is excellent for slow quenches. 



VI. POPULATION INVERSION, DYNAMIC 
CONDUCTIVITY 

Having established the scaling properties of the de- 
fect density in BLG, we turn to the determination of 
the optical response of the excited state resulting from 
the quench, whose momentum distribution is given by 
Eq. ([8]); The occupation number in the upper and lower 
branches of the spectrum is, respectively, f+(p) = P P 
and f-(p) — 1 — P p due to particle- hole symmetry. For 
momenta close to the K point, population inversion oc- 
curs when f+(p) > f-(p), i-e. in the energy range 
2A A <hjj < 2A A \/1 + (7iln2)/(7rAAT), which translates 
in the near adiabatic limit to 

2A A <fiw<2A A + . (15) 

The effect of a small ac electric field can be considered us- 
ing Fermi's golden rule, and the initial dynamic conduc- 
tivity is related to the rate of optical transitions between 
the two states with the same momentum, weighted by the 
probabilities of occupied initial and empty final states, as 

2tt 



T P H = -M^S [ huj - 2^A 2 + e*(p) ) [Up) f + (p)}, 

(16) 

where M p = |w a; (p)eA| is the transition matrix element 
between the higher and lower energy state, where v x (p) = 
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^* + dH / dp x ^i - and A the vector potential. Thence, we 
obtain the dynamic conductivity 



cr(w) = cr 



! exp 



/ 7TT 



-(4A 2 A -(M 5 



(M 5 



4A 2 

J^e(|M-2A A ), 



(17) 



with (T = e 2 /2ft the ac conductivity of BLG-&2 3 .. 

Both absorption and stimulated emission arc taken 
into account, and the negativity of the resulting conduc- 
tivity indicates the dominance of the latter: this indi- 
cates a phase coherent response, which is of course es- 
sential for a laser. In addition, stimulated emission can 
also win against spontaneous emission by increasing the 
intensity of the incoming radiation field. If spontaneous 
emission dominates (luminescence), the resulting radi- 
ation will still be spectrally limited but without phase 
coherence. 

In the frequency range of Eq. 1151 the dynamic conduc- 
tivity is negative due to the population inversion 3 ^ (i.e. 
the energy injected into the system during the quench is 
released) as 



cr(fiw -> 2A A ) w -2cr . 



(18) 



The region of negative conductivity shrinks with increas- 
ing r, without influencing the amplitude of tj(w) precisely 
at the gap edge. This follows from Eq. ([5J), implying 
maximal population inversion at the Dirac point for ar- 
bitrary quench time, i.e. P p= a = 1. For higher frequen- 
cies, <j(lo) is still suppressed with respect to the adiabatic 
optical response. 

The typical lasing frequency lies in the close vicinity 
of Aa (including the THz regime, wavelength of the or- 
der of 10 /im), conveniently tunable by perpendicular 
electric fields^. The relaxation times for intra- and in- 
tcrband processes in MLG are estimated as 1 ps and 1- 
100 ns^i, respectively, which might be further enhanced 
in BLG around half-filling^ 5 .. Thus, the lasing is expected 
to survive for quenching times in the ps-ns range even in 
the presence of the above processes. Repeated quenching 
(like optical pumping) between A and — A is also linked 
to the Kibblc-Zurck theory 3 ^ with similar effects on the 
population inversion. 



The dc conductivity also reveals the effect of the elec- 
tric field quench. In the presence of a clean gap, ex- 
citations, and hence the dc conductivity, are exponen- 
tially suppressed at low temperatures. The excited elec- 
trons/holes in the upper/lower branch resulting from this 
quench, can carry a current that is not activated. 



VII. CONCLUSIONS 

To conclude, by exploiting the tunability of the band 
gap in BLG by a perpendicular electric field, a fi- 
nite rate temporal electric field quench leads to ex- 
cited state production, whose distribution is analyzed 
in terms of Kibble-Zurck scaling, LZ dynamics for non- 
linear quenches and is compared to the full numerical 
solution of the problem with screening corrections, us- 
ing an adiabatic decoupling procedure. The effect of the 
quench is manifested in population inversion, and BLG 
could be used as a coherent sourse of infra-red radiation, 
and possibly as a laser. 

Our results apply to other systems with a quadratic 
band crossing, e.g. for certain nodal superconductors or 
cold atoms on Kagome or checkerboard optical lattices^! 
at appropriate fillings, described by Eq. (TTJ with J = 2 
at low energies. The momentum distribution, Eq. ([5)1 
and the concomitant scaling of the defect density after 
closing and reopening the gap would be direct evidence of 
the quench dynamics. Particularly intriguingly, graphene 
multilayers with appropriate stackings realize higher or- 
der ( J > 2) band crossing a 38 ' 39 . 
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